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The ad vect ion-diffusion equation can be approximated by a one-dimensional diffusion equation in 
Lagrangian coordinates along the directions of compression of fluid elements (the stable manifold). 
This result holds in any number of dimensions, for a velocity field with chaotic trajectories, with 
an error proportional to the square root of the diffusivity. After some time, the one-dimensional 
equation becomes invalid, but by that time a large fraction of the scalar variance has decayed. 
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I. INTRODUCTION 

The importance of transport and mixing processes in 
the physical sciences cannot be overestimated. Applica- 
tions abound in the geophysical sciences such as oceanog- 
raphy, atmospheric science, and geology, as well as in 
astrophysics, chemistry, and general fluid dynamics. Be- 
cause it spans so many disciplines, a comprehensive un- 
derstanding of the detailed manner in which mixing pro- 
ceeds is of paramount importance. 

In most physically relevant applications it is the case 
that the mixing (diffusion) process depends on the stir- 
ring (advection) phase to create large gradients that then 
allow diffusion to act. Without some form of enhanced 
stirring, diffusion is powerless in the face of weak gra- 
dients. The enhancement may be provided by turbu- 
lence in the flow, but that is not necessary: for even 
simple laminar flows, fluid trajectories can display sensi- 
tivity to initial conditions (chaos) that stretches fluid ele- 
ments exponentially — a phenomenon termed chaotic ad- 
vection [1]. This stretching produces the necessary large 
gradients of the advected quantity, leading to chaotic 
mixing. Our focus will be on such smooth, chaotic flows. 

There are two main results in this paper. First, we 
show that the advection-diffusion equation in Lagrangian 
coordinates can be reduced to a one-dimensional diffusion 
equation along the stable manifold of the flow (the direc- 
tion along which fluid elements are compressed). The 
diffusion coefficient of the reduced equation grows expo- 
nentially in time, reflecting the enhancement to diffusion 
due to chaotic trajectories; the reduced equation itself is 
an approximation valid to exponential accuracy in time. 
The enhancement to diffusion is not uniform along the 
stable manifold: the amount of stretching of fluid ele- 
ments is anticorrelated with the curvature of the stable 
manifold [2, 3]. The second result is that this approxi- 
mate, one-dimensional, equation breaks down after a cer- 
tain time because of a buildup of large gradients along 
the neglected direction of stretching. By that time, a sig- 
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nificant fraction of the scalar variance has decayed, so the 
one-dimensional equation is often sufficient for describing 
physical systems. 



II. THE REDUCED ADVECTION-DIFFUSION 
EQUATION 

In the present section we give an account of the man- 
ner in which the advection-diffusion equation expressed 
in Lagrangian coordinates can be collapsed to a one- 
dimensional diffusion equation. The result holds re- 
gardless of dimension (two, three, or more as can arise 
in kinetic equations), and the approximations involved 
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are accurate to order Pe ' , where the Peclet number 
is Pe := Lv/D, with L and v characteristic spatial and 
velocity scales, and D the diffusivity; in many common 
problems, Pe > 10 8 . Though it has been realized for 
some time that the enhanced diffusion proceeds along 
the stable (contracting) direction in both smooth [2, 4] 
and turbulent [5-7] flows, here we make a series of ap- 
proximations and use a geometrical constraint [2, 3, 8] 
to rewrite the advection-diffusion equation as an actual 
one-dimensional diffusion equation along the stable man- 
ifold. This is more than mere formalism, because it al- 
lows the identification of the role of invariant structures 
of the flow in the mixing process (Section III) . Solving 
the one-dimensional equation also allows us to show that 
it must eventually break down (Section IV). 

The advection-diffusion equation for a compressible 
velocity field v(x,t) on an n-dimensional bounded do- 
main is 

^ + «-V^=-V.(pD-V0), (1) 

where the scalar field 4>(x,t) is the concentration of the 
quantity advected, B(jc, t) is the diffusivity tensor, and p 
is the density of the fluid (which satisfies the continuity 
equation) . The advected scalar <fr can represent for exam- 
ple temperature, salinity, or the concentration of a solute 
such as helium. The velocity field v can be obtained by 
solving the Navier-Stokes equation, or it can be a model 
system designed to mimic the qualitative features of a 
realistic flow. The trajectory of a fluid element is given 
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by the solution x(t, to; a) to 
dx 



dt 



= v(x,t), x(t = to, to; a) = a, 



(2) 



so that x(t, to; a) is the coordinate transformation at 
time t from the comoving Lagrangian (a) to stationary 
Eulerian (x) coordinates. For smooth velocity fields, the 
transformation is smooth and invertible. 

By construction, in Lagrangian coordinates the advec- 
tion term drops out, and the advection-diffusion equa- 
tion (1) takes the form 



dt 



i 



V -(p ID ) o-Vo(/>o) 



(3) 



where subscripts denote Lagrangian quantities: 
<t>o{a,t) = (p(x(t,t ;a),t), p (a) = \M\ p(x(t,t ; a),t), 
Vo = 3/ da, and 



»Q = 



-1\T 



(4) 



where the tangent map (or Jacobian) M has compo- 
nents M l q = dx' 1 /da q . The a subscript on the time 
derivative in (3) is a reminder that the derivative is taken 
with a held fixed, as opposed to the derivative in (1) 
which has x fixed. 

The covariance of (1) and (3) (i.e., the form of the 
right-hand side of both equations is the same) suggests 
the definition of a metric 



g := k 1 



(5) 



on the tangent space at x, where the dimensional func- 
tion n(a,t) is chosen such that \g\ = 1. (For constant IB), 
k is also constant.) The metric g is a proper Rieman- 
nian metric, since it is symmetric and positive-definite; 
it induces a metric go on the tangent space at a, 
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(6) 



Because it is symmetric and positive-definite, the met- 
ric go can be written in terms of its real-positive eigen- 
values A 2 and orthonormal eigenvectors e a as 
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<7 = 1 



(7) 



The A a (t,to',(i) are called coefficients of expansion; 
without loss of generality, they are ordered such 
that Ai > A 2 > . . . > A„. For a very wide class of flows, 
including incompressible flows, A„ (the smallest coeffi- 
cient of expansion) is associated with an exponentially 
contracting direction, so that A„ 1. The characteris- 
tic directions e a (t,to;a) converge exponentially to their 
time-asymptotic values e^(to; a) [8, 9]. The direction of 
fastest contraction (a — n) plays a distinguished role in 
our development, so we emphasize its importance by de- 
noting it by the letter 's' (for "stable"). Thus, e„ and A„ 
are written s and A s . Similarly, we will later use the 



letters 'u' (for "unstable") for the direction of fastest 
stretching (cr = 1); for three-dimensional systems, we 
use 'm' (for "median") when referring to direction a = 2, 
which may be stretching, contracting, or neutral. 

Replacing Do by ng^ 1 using (6), and substituting the 
diagonal form (7) of the metric into (3) , we find 



dt 



= E- V ° 



{po^K' 



6<T 6n 



Vo^>) . (8) 



Since k ~ Pe -1 <C 1, for short times the right-hand side 
of (8) can be neglected completely. However, in a chaotic 
flow, at least some of the A CT achieve exponential behavior 
after a moderate time and can overcome the small diffu- 
sivity. Because the inverse of A CT enters (8), the direction 
of fastest contraction A s = A n dominates, and eventually 
we have Pe _1 A s T 2 ~ 0(1). Assuming A s ~ cxp(A s t) this 
occurs roughly at time 



t\ — to — 



2|A S 



log Pe, 



(9) 



after which we can approximate (8) by 



dt 



1 

Po 



Vo • (po«A s - 2 ss • V o 0o) + 0(A r ;! 1 ). (10) 



The relative error we are making if we neglect 
the 0(A~^_ 1 ) terms is of order (A s /A„_i) 2 . Since Pe is as- 
sumed very large, by the time A^ 2 ~ Pe the coefficient of 
expansion A s is well into the exponential regime, and we 
have A s A„_i, if the two coefficients are exponentially 
separated. Thus, for the approximation (10) to hold, it is 
sufficient to require that the two smallest Lyapunov ex- 
ponents differ, which we assume to be the case. In a 3D 
autonomous flow, we have A n _ x = A m ~ 1 because of the 
presence of a zero Lyapunov exponent [10]. Hence, the 
relative error in (10) is of order A 2 - Pe" 1 <C 1. For a 2D 
time-dependent flow, the error is even smaller, being of 
order (A s /A u ) 2 - (Pc _1 / A S) < Pc_1 - We conclude that 
in most physically relevant applications the neglect of the 
higher-order terms in (10) is correct to a very high degree 
of accuracy. However, in Section IV we will show that 
this neglect is not justified for long times. 

Physically, the leading-order term in (10) reflects the 
enhanced diffusion due to the presence of a dominant 
contracting direction: along that direction huge gradi- 
ents are created by the stretching and folding action of 
the flow, leading to an exponential increase in the effi- 
ciency of diffusion. For all practical purposes the other 
directions can be ignored, as long as the smallest Lya- 
punov exponents are nondegenerate, as assumed. 

Since the eigenvector s converges to its asymptotic 
value s°° at a rate (A s /A„_i) [8, 9], we may re- 
place s(t,t ;a) in (10) by s°°(to;a). The error com- 
mitted is proportional to A s ~ Pe -1 / 2 , larger than 
the order Pe -1 error in (10). The characteristic direc- 
tion s°°(to; a) can be integrated to yield the stable man- 
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if old a s (s) (or s-line) through a point a , 
da s (s) 



ds 



= s°°{to;a s (s)), a s (0) = a , (11) 



where s denotes the arc length along the stable manifold. 
The stable manifold is central to our development, and 
we will reformulate (10) so that it represents a diffusion 
equation along this manifold. 

To achieve this goal two more approximations are 
needed. It was shown in Ref. [8] that s • VologA s con- 
verges exponentially to a time- independent function, at a 
rate max(A s , A s /A n _i). This convergence implies that, 
along an s-line given by ct s (s), we can choose an arbi- 
trary reference point a s (0) = a and express the time- 
dependence of A s (t, t a ;a s (s)) as 



A s (t, t ; a) = A s (i ; a) A s (t, t ; o ), 



(12) 



where a and a are on the same s-line and A s is a time- 
independent function. We will discuss the physical sig- 
nificance of the function A s , which we call the parallel 
coefficient of expansion, in Section III. Neglecting the 
time-dependence of A s docs not increase the overall er- 
ror. 

We now rewrite (10) as 



dt 



1 

Po 



V • (p K ss • Vo^o) , 



where we have defined 

k := KA s ~ 2 (t,t ;ao), 



s := 



A^s 00 . 



(13) 



(14) 



In Ref. [8] , it was demonstrated that for a chaotic flow in 
any number of dimensions, s must satisfy the constraint 

div s = A^divos 00 -s°° • V logA s ) =0, (15) 

where 

div s:=| 5o |- 1/2 Vo-(M 1/2 s) (16) 

is the covariant divergence of s [11, p. 49]. Upon using 
the constraint (15) in (13), we obtain 



dt 



where 



-S • V (pK S ■ \7 (j) ) , 

P 



P-=9o 1/2 Po- 



and go(a) is defined analogously to A s in (12), 
\9o(a,t)\ = g (a) \g Q (a ,t)\. 



(17) 



(18) 



(19) 



The decomposition (19) is possible because the derivative 
of the metric determinant along the s-line converges ex- 
ponentially to a time-independent value [8] . The factor go 
in p arises from the metric determinant in the covariant 



divergence (16), and is absent in incompressible flows. 
Thus, after all the approximations have been made, 
the overall error is proportional to max(A s , (A s /A„_i)), 
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which scales as Pe ' . 

Finally, we define a derivative along s, 



9 ~ w 



(20) 



where s is a new parameter along the s-line that differs 
from the arc length s because of the weight A.T 1 in (14). 
Equation (17) then achieves the simple form 



d<t>o 
dt 



1 d 



ds 



(21) 



Equation (21) is a one-dimensional diffusion equation 
along a line defined by a s (s). The effective diffusiv- 
ity k(a s (s),t) diverges exponentially in time, so that 
even for a microscopic physical diffusivity k the gradi- 
ents along s are eventually smoothed. Note that it is 
only in terms of the parameter s that we obtain a one- 
dimensional diffusion equation; in terms of any other 
parameter (other than a uniform rcscaling of s), equa- 
tion (21) contains a "drift" term proportional to V • s°°. 



III. PHYSICAL INTERPRETATION 

We now address the physical meaning of the vari- 
ous constituents of the one-dimensional diffusion equa- 
tion (21). The essence of the simplification in Section II 
lies in expressing the advection-diffusion equation in La- 
grangian coordinates in terms of a single parameter along 
a stable manifold: with an appropriate choice of parame- 
ter, namely s, the equation collapses to a simple diffusion 
equation. 

There are three main ingredients involved in the one- 
dimensional equation (21): (i) The exponential part of 
the coefficient of expansion, A s (t, to; ao); (ii) the stable 
manifold, a s (s); and (iii) The parallel coefficient of ex- 
pansion, A s . We discuss each in turn. 

(i) The exponential part of the coefficient of expan- 
sion, A s (t, to; ao), is responsible for the enhanced dif- 
fusion coefficient k. The variation of the function A s 
across s-lines along unstable directions is very steep (see 
Section IV) and is the cause of the poor convergence of 
Lyapunov exponents [2, 8]. Thus, the diffusion rates on 
two neighboring s-lines can be quite different in mag- 
nitude, though similar in character because A s (i,i ;ao) 
is independent of s and so gives only a time-dependent 
overall scaling of the diffusion coefficient in (21). 

(ii) The s-line is straightforward to compute [2, 3]; it 
has a fine structure and densely fills an crgodic region. 
If we treat the stable manifold as a material line to be 
advected by the flow, then the elements on that line will 
converge exponentially together as t — > oo. Thus, if these 
initial elements carry different concentrations of <j), the 
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gradients of <p will grow exponentially along these tra- 
jectories. It is thus natural for the stable manifold to 
play a central role in the chaotic mixing process. In fact, 
it is the capacity of the stable manifold to fill an entire 
chaotic region that allows for a thorough mixing [12]. 

When evolved backward in time, a small "test parcel" 
(of a scale comparable to the length at which dissipation 
acts) will be stretched and folded and will give a filament 
of fluid that samples the initial data 0o [12]. This fila- 
ment will tend to trace out the stable manifold. This is 
the Lagrangian analogue of the phenomenon of "asymp- 
totic directionality" invoked to discuss the behavior of a 
material line as it is evolved forward in time by a flow [13- 
16]. 

(iii) The parallel coefficient of expansion A s , also 
has deep implications for chaotic mixing. The time- 
dependent part of A s , whose role was explored in (i) 
above, is uniform on an s-line. In contrast, the parallel 
coefficient of expansion A s is constant in time (though 
it depends on t for unsteady flows) but varies along the 
stable manifold. Since the full coefficient of expansion A s 
gives the rate of contraction of an initial infinitesimal el- 
lipsoid, we must regard A s as a deviation from the ref- 
erence rate of contraction given by A s (t, to; ao). Thus, 
regions of anomalously large A s are associated with fluid 
elements that contract slower than other elements on the 
same s-line. 

The size of A s has an impact on diffusion, because 
the enhanced diffusion is due to the contraction of fluid 
elements. The effect is embodied in d/ds: because of (14) 
and (20), we have d/ds = A^ 1 d/ds, where s is the arc 
length along s°°. Hence, in regions of large A s , gradients 
of 4>q are smaller when expressed in the coordinate s; the 
reverse is true for regions of small A s . A large value of A s 
is thus a hindrance to diffusion. 

But what range of values does A s achieve? In Fig. 1 (a) 
we plot logA s against the arc length s on an s-line. The 
cellular flow of Solomon and Gollub [17, 18] is used in the 
computation (with parameter values A=B=u=k= 
1); it is meant to mimic a chain of oscillating convection 
rolls. The parallel coefficient of expansion A s is seen to 
vary over five orders of magnitude. Thus, there exist 
many exceptional fluid elements that have very slow rates 
of contraction. 

The important factor, then, is how frequently these 
large values of A s occur along the s-line. Figure 1(b) 
shows the probability distribution function (PDF) 
of logA s along two s-lincs: the tail of the distribution 
is clearly an exponential, decaying as exp(— 21ogA s ) = 
AjT 2 . This is surprising because for this flow the distri- 
bution of the finite-time Lyapunov exponents themselves 
has a Gaussian tail, as opposed to a "fat" (exponential) 
tail. In fact, we have observed the A~ 2 exponential law 
for the tail of the distribution in several flows, even flows 
such as the random wave flow [4] (also known as the sine 
flow) that have no KAM regions. Thus, the anomalously 



large values of A s cannot only be due to sticking to KAM 
surfaces, though that is certainly a factor. Rather, there 
arc inherent paths in the flow that fluid elements can fol- 
low to avoid stretching, and every "final" fluid element 
contains a mixture of such paths. That the tails of 
the PDF of log A s fall off exponentially rather than as a 
Gaussian suggests that these exceptional events are im- 
portant to the long-time decay of scalar variance, in the 
same manner that the initial alignment of the scalar gra- 
dient along the s-line matters. The statistics of these 
regions of low stretching is related to the distribution of 
curvature of material lines [19, 20]. 



IV. BREAKDOWN OF THE APPROXIMATION 

We now show that the neglect of the higher-order terms 
in (10) is not justified for long times, a fact we believe 
has not been appreciated before. To do that, we compare 
the magnitude of the terms that were neglected in going 
from (8) to (10) and show that they become important af- 
ter a certain time. To simplify the argument, we assume 
the flow is incompressible, and that the diffusivity k is 
constant in space and time, so that k in (14) is constant 
along the s-line and depends only on which particular sta- 
ble manifold we are dealing with [see (i) in Section III]. 
We label the dependence on the particular manifold un- 
der consideration by a parameter u, meant to represent 
displacements along an unstable (stretching) direction u 
of the flow. For an initial condition 4>o(s, 0; u) — exp(ifcs) 
(we set the initial time to to zero), the one-dimensional 
diffusion equation (21) has solution 

cj> (s,t;u) = cxp (iks - k 2 tX{t;u)) , (22) 

where 



X{t;u) := - / k{t';u)At' (23) 
t Jo 

is the time-average of k(t;u). From Eq. (14), k is pro- 
portional to AjT 2 ~ exp(2|A s |t), so (22) gives a superex- 
ponential decay of (f> a . 

Consider now the variation of <j)o(s, t; u) along the un- 
stable direction u, 



d u <t>v = -k 2 td u X(t;u) 4> , 



(24) 



where d u := u • Vq. The second derivative of 0o(s, t; u) is 



(k 2 t d u X(t; u)) 2 ~k 2 t d 2 u X{t- u) 



(25) 



Keeping only the first term in (25) (the second term does 
not dominate), we take the ratio of the a = 1 (unsta- 
ble direction) term to the a = 3 (stable direction) term 
in (8), 



s := 



| ac A" 2 {k 2 td u X) 2 fa\ 
Ik Ag~ 2 fc 2 0o| 



= \A- 1 A s kt(d u X)\ 2 ; (26) 
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(a) 



(b) 



FIG. 1: (a) log A s as a function of arc length s along a typical s-line. (b) Normalized PDF of A s on two s-lines of equal length. 
For comparison, the dashed line denotes A~ 2 . The PDFs on the two s-lines appear very similar. 



The ratio e must be small to justify the approximation 
leading from (8) to (10). Because k grows exponentially, 
we have X ~ «;A~ 2 /2|A s |i, so that 



d u X 



_ q A" 2 - — 
2|A s |i " As |A s |i 



A^^logA" 1 . 



(27) 



Since A s is in the exponential regime, we have that [8] 

d„logA s = cA u , (28) 

where c is a function of space and has only algebraic 
behavior in time: it neither grows nor decays exponen- 
tially in time. The exponential growth in (28) is a re- 
flection of the exponential separation of fluid trajectories 
in the chaotic flow. Equation (28) applies for almost 
all initial conditions in essentially all smooth flows and 
maps (it is sufficient that the velocity field be thricc- 
differentiable [8]), except in highly exceptional cases such 
as uniformly stretching systems (e.g., Arnold's cat map, 
or flow around a hyperbolic point), for which e = 0; the 
superexponential solution (22) then persists forever. 
After inserting (27)-(28) in (26), we obtain finally 



\kc\ a k A a 



Pe A; 



(29) 



Assuming that A s 1 ~ exp(|A s |i), we conclude that e be- 
comes 0(1) at time 



t 2 ~ -i- logPc = 2ii 



(30) 



where t\ is the time where diffusion sets in, defined by (9). 

In three dimensions, we must take into account 
the cr = 2 (the 'm' or 'median' direction) term in (8). 



Its role depends on the type of direction it represents. 
If it is a stretching or neutral direction, then the esti- 
mate presented above holds, and diffusion along the a = 
2 direction becomes important at roughly the same 
time t 2 given by (30). If a = 2 is a contracting di- 
rection, we must modify the asymptotic behavior (28) 
to <9 m logA s = c(A s /A m ), as described in [8]. This leads 
to a time t' 2 ■= (l/2|A m |) logPe, after which the a = 2 
term in (8) can no longer be neglected. Depending on 
the magnitude of |A m |, t 2 may come before or after t 2 . 
The breakdown of (21) thus occurs at min^,^), which 
we will simply call t 2 . In two dimensions, we ignore t' 2 
because there is no median direction. 

The overall picture is thus as follows. For time t <t\, 
0o is approximatively constant because the exponential 
growth of gradients has not yet overcome the small diffu- 
sivity. For t\ < t <t 2 , Eq. (21) applies and the scalar 
variance decays superexpontially. For t > t 2 , the ne- 
glect of the (7 = 1 (and a = 2 in 3D) terms in (8) is 
no longer justified, and the one-dimensional diffusion 
equation (21) must be abandoned: to characterize the 
advection-diffusion process, the full Eq. (8) must be 
solved. 

The reason that the diffusion along u cannot be ne- 
glected after t 2 is that the tiny, exponentially-decaying 
effective diffusivity associated with u allows the creation 
of extremely large gradients of (f> along that direction. 
These gradients are a consequence of the superexponen- 
tial solution (22). Eventually, the gradients overcome the 
exponentially-decaying effective diffusivity. In essence, 
the superexponential solution contains the seeds of its 
own destruction: the superexponentiality comes from the 
factor X in (22), but that factor exhibits extremely rapid 
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variation along the stretching direction, with no diffusiv- 
ity to oppose it (until t ~ t 2 ). 

The breakdown of the one-dimensional equation helps 
to understand a discrepancy between the superexponen- 
tial decay predicted by that equation and the observed 
exponential decay of the scalar variance in typical chaotic 
and turbulent flows. In Rcf. [4] a superexponential de- 
cay was derived for individual fluid trajectories, and an 
exponential decay was recovered by averaging over the 
initial angle between the scalar gradient and the stable 
direction. Fereday et al. [21] pointed out that the decay 
of the variance in the (nonuniform) baker's map is ex- 
ponential for long times, even though in their case the 
angle between the scalar gradient and the stable direc- 
tion is always zero, so that no averaging can be done. 
Our demonstration of the termination of the superexpo- 
nential regime after some time shows that there is hope 
of restoring the more appropriate, exponential decay of 
the variance, but we have yet to show this because a tech- 
nique of solution for the Eq. (8) must be found. For the 
Kraichnan-Kazantsev model of homogeneous, isotropic 



turbulence (delta-correlated in time), the decay of vari- 
ance is known to be exponential [5-7], though the pre- 
dicted rate is too rapid when compared with nonhomogc- 
neous systems [21] which can have slowly-decaying eigen- 
functions. 

Even though it ceases to be valid after some time, the 
one-dimensional equation is still relevant, because during 
the time t\ < t < t<i a very large fraction of the variance 
can get dissipated, owing to the very fast decay of the 
superexponential. 
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